//
// perc2 - calculation of the fractal dimension of correlated
// bond percolation cluster hulls
//
// Copyright (C) 2009, 2010 Indrek Mandre <indrek(at)mare.ee>
// http://www.mare.ee/indrek/perc2/, http://code.google.com/p/perc2/
// 
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
// 

#include <sstream>

#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <sys/mman.h>

#include <string.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include <sys/types.h>
#include <sys/stat.h>

#include "perc2.hpp"

#ifdef GMP_VALUE
  class _gmp_value_init
  {
  public:
    _gmp_value_init ()
    {
      mpf_set_default_prec (GMP_PREC);
    }
  };
  static _gmp_value_init __gmp_value_init_obj;
#endif

double EXPONENT = 0.75;
double BASEPOW = DISTPOWER;
int LAYERSTART = -1;
int LAYEREND = LAYERCOUNT;
RNDF rndf;

value_t calc_cell () __PERC_INLINE;

value_t calc_cell ()
{
  value_t px = rndf.rnd1();
  value_t py = rndf.rnd1();
  value_t add = 0;
  //double cx, cy, d2;
  for ( int i = -1; i <= 1; i++ )
    {
      for ( int j = -1; j <= 1; j++ )
        {
          value_t bcx = 0.5 + i;
          value_t bcy = 0.5 + j;
          value_t dist2 = (px - bcx) * (px - bcx) + (py - bcy) * (py - bcy);
          value_t hf = dist2hf (dist2);
          if ( hf != 0 )
              add += rndf.rnd() * hf;
        }
    }
  return add;
}

MemPool::MemPool (size_t size, size_t grow) : size(size), grow(grow), chunks(0), last(0)
{
  extend();
  reset();
}

MemPool::~MemPool ()
{
  while ( chunks )
    {
      char *mem = (char *)chunks;
      chunks = chunks->next;
      delete [] mem;
    }
}

void MemPool::extend ()
{
  link_t *c = (link_t *)new char [MEMPOOL_LINK_SIZE + grow * size];
  c->next = 0;
  if ( last )
      last->next = c;
  else
      chunks = c;
  last = c;
}

void MemPool::reset ()
{
  current = chunks;
  begin = ((char *)current) + MEMPOOL_LINK_SIZE;
  end = begin + size * grow;
}

MemPool * Chunk::pool = new MemPool (sizeof (Chunk), 4096);

Field::Field (bool add_baseline) : add_baseline (add_baseline)
#ifndef UNCORRELATED
  , db0(1 + FIELD_SIZE / CHUNK_SIZE)
#endif
{
#if defined(FORCE_HEIGHT) || defined(UNCORRELATED)
  layers.push_back (new Layer (0));
#else
  int lowest_layer = LAYEREND;
  while (1)
    {
      lowest_layer--;
      if ( lowest_layer == LAYERSTART - 1 )
          break;
      if ( pow(BASEPOW, lowest_layer) * 3 <= 1 )
          break;
    }

  for ( int l = LAYEREND; l > lowest_layer; l-- )
    {
      for ( size_t j = 0; j < CUMULATION; j++ )
          layers.push_back (new Layer (l));
    }

  baserng.init (lowest_layer - (lowest_layer - LAYERSTART), lowest_layer);
#endif
  reset ();
}

value_t Field::value (double xc, double yc)
{
#ifdef UNCORRELATED
  return layers[0]->value (xc, yc);
#endif
  double dx = FIELD_SIZE / 2 - xc;
  double dy = FIELD_SIZE / 2 - yc;
  value_t ret = 0;
  for ( size_t i = 0; i < layers.size(); i++ )
      ret += layers[i]->value (dx, dy);
#ifndef UNCORRELATED
  ret += db0.get_by_layer (xc, yc, baserng);
  ret += baselinernd;
#endif
  return ret;
}

void Field::printsurf (const char *path, int width, int height, bool boolean)
{
  FILE *fp = fopen (path, "w");
  for ( int y = 0; y <= height; y++ )
    {
      for ( int x = 0; x <= width; x++ )
        {
          //double v = value (midpos() - width/2 + x, midpos() - height/2 + y);
          double v = v2d(value (x, y));
          if ( boolean )
              v = (v < 0) ? -1 : 1;
          fprintf (fp, "%.3g ", v);
        }
      fprintf (fp, "\n");
    }
  fclose (fp);
}

class fbitmap
{
public:
  fbitmap (size_t dim = FIELD_SIZE)
  {
    ysize = dim + 2;
    xsize = 2 + (dim >> 3);

    data = new char [xsize * ysize];
    memset (data, 0, xsize * ysize);
  }

  ~fbitmap ()
  {
    delete [] data;
  }

  inline void set (int x, int y) __PERC_INLINE
  {
    at(x,y) |= 1 << (x & 0x7);
  }

  inline bool get (int x, int y) const __PERC_INLINE
  {
    return at(x, y) & (1 << (x & 0x7));
  }

  void clear ()
  {
    memset (data, 0, xsize * ysize);
  }

private:
  size_t ysize, xsize;
  char *data;

  fbitmap (const fbitmap&);
  const fbitmap& operator= (const fbitmap&);

  char& at (int x, int y) __PERC_INLINE
  {
    return data[y * xsize + (x >> 3)];
  }

  const char& at (int x, int y) const __PERC_INLINE
  {
    return data[y * xsize + (x >> 3)];
  }
};

static inline size_t is_single (const fbitmap& valmap, int x, int y)
{
  if ( Field::upw (x, y) )
    {
      if ( (!valmap.get (x, y + 1) &&
          !valmap.get (x + 1, y) &&
          !valmap.get (x + 1, y + 1)) ||
          (!valmap.get (x, y - 1) &&
          !valmap.get (x - 1, y) &&
          !valmap.get (x - 1, y - 1)) )
          return 1;
    }
  else
    {
      if ( (!valmap.get (x, y + 1) &&
          !valmap.get (x + 1, y) &&
          !valmap.get (x + 1, y - 1)) ||
          (!valmap.get (x, y - 1) &&
          !valmap.get (x - 1, y) &&
          !valmap.get (x - 1, y + 1)) )
          return 1;
    }
  return 0;
}

static inline size_t is_full (const fbitmap& valmap, int x, int y)
{
  if ( Field::upw (x, y) )
    {
      if (
          valmap.get (x, y + 1) &&
          valmap.get (x + 1, y) &&
          valmap.get (x + 1, y + 1) &&
          valmap.get (x, y - 1) &&
          valmap.get (x - 1, y) &&
          valmap.get (x - 1, y - 1) )
          return 1;
    }
  else
    {
      if (
          valmap.get (x, y + 1) &&
          valmap.get (x + 1, y) &&
          valmap.get (x + 1, y - 1) &&
          valmap.get (x, y - 1) &&
          valmap.get (x - 1, y) &&
          valmap.get (x - 1, y + 1) )
          return 1;
    }
  return 0;
}

static inline size_t is_line (const fbitmap& valmap, int x, int y)
{
  if ( Field::upw (x, y) )
    {
      if (
          !valmap.get (x, y + 1) &&
          !valmap.get (x + 1, y) &&
          valmap.get (x + 1, y + 1) &&
          !valmap.get (x, y - 1) &&
          !valmap.get (x - 1, y) &&
          valmap.get (x - 1, y - 1) )
          return 1;
    }
  else
    {
      if (
          !valmap.get (x, y + 1) &&
          !valmap.get (x + 1, y) &&
          valmap.get (x + 1, y - 1) &&
          !valmap.get (x, y - 1) &&
          !valmap.get (x - 1, y) &&
          valmap.get (x - 1, y + 1) )
          return 1;
    }
  return 0;
}

// Calculate the clusters' fractal dimension
void Field::testdim ()
{
#if 1
  int midx = FIELD_SIZE / 2 + 1;
  int midy = FIELD_SIZE / 2 + 1;

  std::vector<aggregator<5> > res;
  res.resize (FIELD_SIZE);

#define CCC 10000
  for ( size_t i = 0; i < CCC; i++ )
    {
      while (1)
        {
          reset();
          if ( valueb(midx, midy) )
              break;
        }
      fbitmap *scanmap = new fbitmap ();
      fbitmap *valmap = new fbitmap ();

      std::vector<std::pair<int,int> > lookinto;
      std::vector<std::pair<int,int> > buildup;

      lookinto.push_back (std::pair<int,int>(midx, midy));

      size_t cluster_size = 0;
      int have_lu, have_lb, have_ru, have_rb, have_tl, have_tr, have_bl, have_br;

      have_lu = have_lb = have_ru = have_rb = have_tl = have_tr = have_bl = have_br = 0;

      while(1)
        {
          buildup.clear();
          bool added = false;
          for ( size_t i = 0; i < lookinto.size(); i++ )
            {
              int x = lookinto[i].first;
              int y = lookinto[i].second;
              if ( !scanmap->get (x, y) )
                {
                  scanmap->set (x, y);

                  if ( x <= 0 )
                    {
                      if ( y >= FIELD_SIZE / 2 )
                          have_lu++;
                      else
                          have_lb++;
                      continue;
                    }
                  if ( x >= FIELD_SIZE )
                    {
                      if ( y >= FIELD_SIZE / 2 )
                          have_ru++;
                      else
                          have_rb++;
                      continue;
                    }
                  if ( y <= 0 )
                    {
                      if ( x >= FIELD_SIZE / 2 )
                          have_tr++;
                      else
                          have_tl++;
                      continue;
                    }
                  if ( y >= FIELD_SIZE )
                    {
                      if ( x >= FIELD_SIZE / 2 )
                          have_br++;
                      else
                          have_bl++;
                      continue;
                    }

                  if ( valueb (x, y) )
                    {
                      valmap->set (x, y);
                      added = true;
                      cluster_size++;
                      if ( upw (x, y) )
                        {
                          buildup.push_back (std::pair<int,int>(x - 1, y));
                          buildup.push_back (std::pair<int,int>(x + 1, y));
                          buildup.push_back (std::pair<int,int>(x, y + 1));
                          buildup.push_back (std::pair<int,int>(x, y - 1));
                          buildup.push_back (std::pair<int,int>(x - 1, y - 1));
                          buildup.push_back (std::pair<int,int>(x + 1, y + 1));
                        }
                      else
                        {
                          buildup.push_back (std::pair<int,int>(x - 1, y));
                          buildup.push_back (std::pair<int,int>(x + 1, y));
                          buildup.push_back (std::pair<int,int>(x, y + 1));
                          buildup.push_back (std::pair<int,int>(x, y - 1));
                          buildup.push_back (std::pair<int,int>(x - 1, y + 1));
                          buildup.push_back (std::pair<int,int>(x + 1, y - 1));
                        }
                    }
                }
            }
          if ( !added )
              break;
          lookinto.swap (buildup);
        }

      //if ( !(have_l > FIELD_SIZE / 16 && have_r > FIELD_SIZE / 16 && have_t > FIELD_SIZE / 16 && have_b > FIELD_SIZE / 16) )
      #define TC 1
      printf ("Testing, cluster_size=%zu, parts=%d, %d, %d, %d, %d, %d, %d, %d\n", cluster_size, have_lu, have_lb, have_ru, have_rb, have_tl, have_tr, have_bl, have_br);
      if ( !(have_lu > TC && have_lb > TC &&
          have_ru > TC && have_rb > TC &&
          have_tl > TC && have_tr > TC &&
          have_bl > TC && have_br > TC) )
        {
          printf ("Fail!\n");
          i--;
          delete scanmap;
          delete valmap;
          continue;
        }

      printf ("Recording, i=%zu\n", i);

      size_t c = 0;
      int s = 0;
      size_t singles = 0, fulls = 0, lines = 0;
      while ( s < 512 )
        {
          for ( int x = midx - s; x <= midx + s; x++ )
            {
              int y1 = midy - s;
              int y2 = midy + s;
              if ( valmap->get (x, y1) )
                  c++;
              if ( valmap->get (x, y2) )
                  c++;

              singles += is_single(*valmap, x, y1);
              singles += is_single(*valmap, x, y2);
              fulls += is_full(*valmap, x, y1);
              fulls += is_full(*valmap, x, y2);
              lines += is_line(*valmap, x, y1);
              lines += is_line(*valmap, x, y2);
            }

          for ( int y = midy - s + 1; y <= midy + s - 1; y++ )
            {
              int x1 = midx - s;
              int x2 = midx + s;
              if ( valmap->get (x1, y) )
                  c++;
              if ( valmap->get (x2, y) )
                  c++;
              singles += is_single(*valmap, x1, y);
              singles += is_single(*valmap, x2, y);
              fulls += is_full(*valmap, x1, y);
              fulls += is_full(*valmap, x2, y);
              lines += is_line(*valmap, x1, y);
              lines += is_line(*valmap, x2, y);
            }
          res[2*s+1].add (c, singles, fulls, lines, 0);
          s++;
        }

      delete scanmap;
      delete valmap;

      FILE *fp = fopen ("clscale.csv", "w");
      assert (fp);
      for ( size_t j = 0; j < FIELD_SIZE; j++ )
        {
          if ( !res[j].empty() )
              fprintf (fp, "%zu,%s\n", j, res[j].sprint());
              //fprintf (fp, "%zu %g\n", j, (double)res[j]/(i + 1));
        }
      fclose (fp);
    }
#endif

#if 0
  int midx = FIELD_SIZE / 2 + 1;
  int midy = FIELD_SIZE / 2 + 1;
  int ITC = 16777216;
  for ( int S = 16; S <= 2048; S *= 2 )
    {
      size_t cumsize = 0;
      for ( int i = 0; i < ITC; i++ )
        {
          while (1)
            {
              reset();
              if ( valueb(midx, midy) )
                  break;
            }
          fbitmap *scanmap = new fbitmap ();

          std::vector<std::pair<int,int> > lookinto;
          std::vector<std::pair<int,int> > buildup;

          lookinto.push_back (std::pair<int,int>(midx, midy));

          size_t cluster_size = 0;

          while(1)
            {
              buildup.clear();
              bool added = false;
              for ( size_t i = 0; i < lookinto.size(); i++ )
                {
                  int x = lookinto[i].first;
                  int y = lookinto[i].second;
                  if ( abs (x - midx) > S || abs (y - midy) > S )
                      continue;
                  if ( !scanmap->get (x, y) )
                    {
                      scanmap->set (x, y);
                      if ( valueb (x, y) )
                        {
                          added = true;
                          cluster_size++;
                          buildup.push_back (std::pair<int,int>(x - 1, y));
                          buildup.push_back (std::pair<int,int>(x + 1, y));
                          buildup.push_back (std::pair<int,int>(x, y + 1));
                          buildup.push_back (std::pair<int,int>(x, y - 1));
                          buildup.push_back (std::pair<int,int>(x - 1, y - 1));
                          buildup.push_back (std::pair<int,int>(x + 1, y + 1));
                        }
                    }
                }
              if ( !added )
                  break;
              lookinto.swap (buildup);
            }

          //printf ("Cluster size: for box %d: %zu\n", S, cluster_size);

          delete scanmap;
          cumsize += cluster_size;
        }
      fprintf (fp, "%d %g\n", 2*S+1, (double)cumsize / ITC);
      fflush (fp);
      printf ("%d %g\n", 2*S+1, (double)cumsize / ITC);
      ITC /= 2;
    }
#endif

#if 0
  fbitmap *scanmap = new fbitmap ();
  fbitmap *valmap = new fbitmap ();
  int midx = FIELD_SIZE / 2;
  int midy = FIELD_SIZE / 2;
  while ( !valueb(midx, midy) )
      reset();

  std::vector<std::pair<int,int> > lookinto;
  std::vector<std::pair<int,int> > buildup;

  lookinto.push_back (std::pair<int,int>(midx, midy));

  size_t cluster_size = 0;

  while(1)
    {
      buildup.clear();
      bool added = false;
      for ( size_t i = 0; i < lookinto.size(); i++ )
        {
          int x = lookinto[i].first;
          int y = lookinto[i].second;
          if ( x < 0 || y < 0 || x > FIELD_SIZE || y > FIELD_SIZE )
              continue;
          if ( !scanmap->get (x, y) )
            {
              scanmap->set (x, y);
              if ( valueb (x, y) )
                {
                  valmap->set (x, y);
                  added = true;
                  cluster_size++;
                  buildup.push_back (std::pair<int,int>(x - 1, y));
                  buildup.push_back (std::pair<int,int>(x + 1, y));
                  buildup.push_back (std::pair<int,int>(x, y + 1));
                  buildup.push_back (std::pair<int,int>(x, y - 1));
                  buildup.push_back (std::pair<int,int>(x - 1, y - 1));
                  buildup.push_back (std::pair<int,int>(x + 1, y + 1));
                }
            }
        }
      if ( !added )
          break;
      lookinto.swap (buildup);
    }

  printf ("Cluster size: %zu\n", cluster_size);

  FILE *fp = fopen ("clscale.csv", "w");
  assert (fp);

  size_t c = 1;
  int s = 1;
  while ( s < FIELD_SIZE / 2 - 1 )
    {
      for ( int x = FIELD_SIZE / 2 - s; x <= FIELD_SIZE / 2 + s; x++ )
        {
          if ( valmap->get (x, FIELD_SIZE / 2 - s) )
              c++;
          if ( valmap->get (x, FIELD_SIZE / 2 + s) )
              c++;
        }

      for ( int y = FIELD_SIZE / 2 - s + 1; y <= FIELD_SIZE / 2 + s - 1; y++ )
        {
          if ( valmap->get (FIELD_SIZE / 2 - s, y) )
              c++;
          if ( valmap->get (FIELD_SIZE / 2 + s, y) )
              c++;
        }

      if ( s == 2 || s == 4 || s == 8 || s == 16 || s == 32 || s == 64 || s == 128 || s == 256 || s == 512 )
          fprintf (fp, "%d %zu\n", 2 * s + 1, c);
      s++;
    }

  delete scanmap;
  delete valmap;

  fclose (fp);
#endif
}

void Field::reset ()
{
  //printf ("reset!\n");
  for ( size_t i = 0; i < layers.size(); i++ )
      layers[i]->reset();
#ifndef UNCORRELATED
  db0.reset();
#endif
  Chunk::reset();
  //printf ("baseline: %g\n", sqrt((BASELINECORRECTION * ALPHA * pow(2, -EXPONENT * (LAYEREND + 1)))));
  //exit(0);

#ifndef UNCORRELATED
  if ( add_baseline )
    {
      baseline = BASELINECORRECTION * ALPHA * vpow(BASEPOW, -EXPONENT * (LAYEREND + 1));

      baselinernd = 0;
      for ( size_t i = 0; i < CUMULATION; i++ )
        {
    #if 1
          baselinernd += vsqrt(baseline) * calc_cell();
    #else
          for ( size_t j = LAYEREND + 1; j <= 800; j++ )
            {
              value_t h = vsqrt ( ALPHA * (1 - vpow (BASEPOW, -EXPONENT))) * vpow (BASEPOW, -EXPONENT * (j / 2.0));
              baselinernd += h * calc_cell();
            }
    #endif
        }
    }
  else
    {
      baseline = 0;
      baselinernd = 0;
    }
#endif

  sealevel = 0;
}

void Field::reset_offsets ()
{
  for ( size_t i = 0; i < layers.size(); i++ )
      layers[i]->reset_offset();
}

class Bitmap
{
public:
  Bitmap (int size = 32)
  {
    side = -1;
    cnt = 0xff;
    data = 0;
    reset(size);
  }

  ~Bitmap ()
  {
    free (data);
  }

  void reset (int size)
  {
    if ( side != 2 * size + 4 )
      {
        side = 2 * size + 4;
        hside = side / 2;
        if ( data )
            free (data);
        data = (unsigned char *)calloc (side * side, 1);
        cnt = 1;
        return;
      }

    cnt++;

    if ( !cnt )
      {
        memset (data, 0, side * side);
        cnt = 1;
      }
  }

  inline bool set (int x, int y) __PERC_INLINE
  {
    x += hside;
    y += hside;
    if ( data[y * side + x] == cnt )
        return false;
    data[y * side + x] = cnt;
    return true;
  }

  inline bool get (int x, int y) __PERC_INLINE
  {
    x += hside;
    y += hside;
    return data[y * side + x] == cnt;
  }
private:
  unsigned char *data;
  unsigned char cnt;
  int side, hside;
};

#if 0
class Bitmap
{
  #define SUBSET_LEVEL 6
  #define SUBSET_SIDE (1<<SUBSET_LEVEL)
  #define SUBSET_BITMAP 0x3f
  struct subset
  {
    unsigned char data[SUBSET_SIDE][(SUBSET_SIDE >> 3) + 1];

    inline subset () __PERC_INLINE
    {
      memset (data, 0, sizeof(data));
    }

    inline bool set (int x, int y) __PERC_INLINE
    {
      int bi = x & 0x7;
      unsigned char& d = data[y][x >> 3];
      if ( d & (1 << bi) )
          return false;
      d |= (1 << bi);
      return true;
    }

    inline bool get (int x, int y) __PERC_INLINE
    {
      int bi = x & 0x7;
      unsigned char& d = data[y][x >> 3];
      if ( d & (1 << bi) )
          return true;
      return false;
    }

    static inline void * operator new (size_t size) __PERC_INLINE
    {
      return pool->get();
    }

    static inline void operator delete (void *mem, size_t)
    {
      if ( mem )
          pool->release (mem);
    }

    static MemPool *pool;
  };

public:
  Bitmap (int size)
  {
    reset (size);
  }

  void reset (int size)
  {
    offset = size + 2;
    sidelen = offset * 2;
    ssc = (sidelen / SUBSET_SIDE) + 1;
    data.clear();
    data.resize (ssc * ssc, 0);
    subset::pool->reset();
  }

  inline bool set (int x, int y) __PERC_INLINE
  {
    x += offset;
    y += offset;
    size_t coord = ssc * (y >> SUBSET_LEVEL) + (x >> SUBSET_LEVEL);
    if ( !data[coord] )
        data[coord] = new subset();
    return data[coord]->set (x & SUBSET_BITMAP, y & SUBSET_BITMAP);
  }

  inline bool get (int x, int y) __PERC_INLINE
  {
    x += offset;
    y += offset;
    size_t coord = ssc * (y >> SUBSET_LEVEL) + (x >> SUBSET_LEVEL);
    if ( !data[coord] )
        data[coord] = new subset();
    return data[coord]->get (x & SUBSET_BITMAP, y & SUBSET_BITMAP);
  }

private:
  std::vector<subset *> data;
  int offset;
  int sidelen;
  int ssc;
};

MemPool * Bitmap::subset::pool = new MemPool (sizeof (Bitmap::subset), 256);
#endif

static inline unsigned char makeway(int b1, int b2, int b3, int b4, int b5, int b6, int b7, int b8) __PERC_INLINE;
static inline unsigned char makeway(int b1, int b2, int b3, int b4, int b5, int b6, int b7, int b8)
{
  return
      (b1 << 7) +
      (b2 << 6) +
      (b3 << 5) +
      (b4 << 4) +
      (b5 << 3) +
      (b6 << 2) +
      (b7 << 1) +
      (b8 << 0);
}

int Walker::prewalk (int limit)
{
  bool ret = false;
  coord_t cur = start;

  static Bitmap bmap(32);
  bmap.reset (limit + 1);

  while ( 1 )
    {
      int ox = start.pos.x - cur.pos.x;
      int oy = start.pos.y - cur.pos.y;

      if ( abs (ox) == limit || abs (oy) == limit )
        {
          ret = true;
          break;
        }
          
      if ( field->valueb (cur.pos.x, cur.pos.y) )
        {
          bmap.set (cur.pos.x - start.pos.x, cur.pos.y - start.pos.y);
          cur.hamster();
        }
      else
        {
          cur.punch();
          cur.hamster();
        }
      if ( cur == start )
          break;
    }

  int ox = start.pos.x - cur.pos.x;
  int oy = start.pos.y - cur.pos.y;
  return std::max (abs (ox), abs (oy));
}

bool Walker::walk (int limit)
{
  //int limit2 = limit * limit;
  bool ret = false;
  coord_t cur = start;

#ifdef PRINTOUT
  path_t path;
  grid_t grid;
  path.push_back (cur);
#endif

  hull = 0;
  bonds = 0;
  ends = 0;
  lines = 0;
  sides = 0;
  fjords = 0;
  corners = 0;
  int turns = 0;
  int lc = 0;
  bool was_left = true;
  static Bitmap bmap(32);
  bmap.reset (limit + 1);
  unsigned char pathway = 0;

  while ( 1 )
    {
      int ox = start.pos.x - cur.pos.x;
      int oy = start.pos.y - cur.pos.y;
      int x = ox;
      int y = oy;

      if ( abs (x) == limit || abs (y) == limit )
        {
          ret = true;
          if ( lc >= 4 )
              lines += lc - 3;
          break;
        }
          
#if 0
      int dist2 = (start.pos.x - cur.pos.x) * (start.pos.x - cur.pos.x) +
          (start.pos.y - cur.pos.y) * (start.pos.y - cur.pos.y);
      if ( dist2 >= limit2 )
        {
          len = sqrt (dist2);
          ret = true;
          if ( lc >= 4 )
              lines += lc - 3;
          break;
        }
#endif

      if ( field->valueb (cur.pos.x, cur.pos.y) )
        {
          pathway = (pathway << 1) + 1;
          if ( bmap.set (cur.pos.x - start.pos.x, cur.pos.y - start.pos.y) )
              bonds++;
          else
              sides++;
          // bond, we are doing a left turn
          if ( was_left )
            {
              if ( lc >= 4 )
                  lines += lc - 3;
              lc = 1;
            }
          else
            {
              lc++;
              was_left = true;
            }

          //bonds++;
          turns = 1;
#ifdef PRINTOUT
          grid.push_back (cur.pos);
#endif
          cur.hamster();
#ifdef PRINTOUT
          path.push_back (cur);
#endif
        }
      else
        {
          pathway = (pathway << 1) + 0;
          // no bond, we are doing a right turn
          if ( was_left )
            {
              lc++;
              was_left = false;
            }
          else
            {
              if ( lc >= 4 )
                  lines += lc - 3;
              lc = 1;
            }

          cur.punch();
#ifdef PRINTOUT
          path.push_back (cur);
#endif
          cur.hamster();
#ifdef PRINTOUT
          path.push_back (cur);
#endif
          turns++;
          if ( turns == 4 )
              ends++;
        }
      if ( cur == start )
          break;
      if ( pathway == makeway (0,1,0,1,1,1,0,1) )
          fjords++;
      if ( pathway == makeway (0,1,0,1,1,0,1,0) )
          corners++;
      hull++;
    }

#ifdef PRINTOUT
  printf ("Printout: printing now, having taken %zu hull steps\n", hull);
  print_path ("printout_path.csv", path);
  print_grid ("printout_grid.csv", grid);
  print_surf (grid);
#endif
  return ret;
}

bool Walker::walk_perim2 (int limit)
{
  bool ret = false;

#ifdef PRINTOUT
  path_t path, path1, path2;
  grid_t grid;
  std::vector<pos_t> perimpath;
#endif

  perim = bonds = ends = lines = holes = fjords = corners = 0;
#if 1
  int turns = 0;
  int lc = 0;
  bool zigzag = true;
  bool corner = false;
#endif
  static Bitmap bmap(32);
  bmap.reset (limit + 1);
  int hlimit = limit / 2;

  coord_t cur = start;

#ifdef PRINTOUT
  path1.push_back (cur);
#endif

  while ( 1 )
    {
      if ( abs (start.pos.x - cur.pos.x) == hlimit || abs (start.pos.y - cur.pos.y) == hlimit )
          break;
          
      if ( field->valueb (cur.pos.x, cur.pos.y) )
        {
#ifdef PRINTOUT
          grid.push_back (cur.pos);
#endif
          bmap.set (cur.pos.x - start.pos.x, cur.pos.y - start.pos.y);
          cur.hamster();
#ifdef PRINTOUT
          path1.push_back (cur);
#endif
        }
      else
        {
          cur.punch();
#ifdef PRINTOUT
          path1.push_back (cur);
#endif
          cur.hamster();
#ifdef PRINTOUT
          path1.push_back (cur);
#endif
        }
      if ( cur == start )
          return false;
    }

  coord_t end = cur;

  coord_t end1 = cur;

  cur = start;

#ifdef PRINTOUT
  path2.push_back (cur);
#endif

  coord_t save;

  while ( 1 )
    {
      if ( abs (start.pos.x - cur.pos.x) == hlimit || abs (start.pos.y - cur.pos.y) == hlimit )
          break;
          
      if ( field->valueb (cur.pos.x, cur.pos.y) )
        {
          save = cur;
#ifdef PRINTOUT
          grid.push_back (cur.pos);
#endif
          bmap.set (cur.pos.x - start.pos.x, cur.pos.y - start.pos.y);
          cur.rhamster();
#ifdef PRINTOUT
          path2.push_back (cur);
#endif
        }
      else
        {
          cur.punch();
#ifdef PRINTOUT
          path2.push_back (cur);
#endif
          cur.rhamster();
#ifdef PRINTOUT
          path2.push_back (cur);
#endif
        }
      if ( cur == start )
          return false;
    }

  if ( abs(end1.pos.x - cur.pos.x) < limit &&
      abs(end1.pos.y - cur.pos.y) < limit )
      return false;

  coord_t end2 = save;

  cur = end2;

#ifdef PRINTOUT
  path.push_back (cur);
#endif

  while (1)
    {
#if 0
      int ox = end2.pos.x - cur.pos.x;
      int oy = end2.pos.y - cur.pos.y;
      int x = ox;
      int y = oy;

      //if ( abs (x) == limit || abs (y) == limit )
#endif
      if ( cur == end1 )
        {
          ret = true;
          if ( lc >= 6 )
              lines += lc - 5;
          break;
        }

      int rx = cur.pos.x - start.pos.x;
      int ry = cur.pos.y - start.pos.y;

      bool bounce = bmap.get (rx, ry);

      if ( bounce )
          bonds++;

      if ( !bounce )
        {
          if ( (cur.pos.x ^ cur.pos.y) & 1 )
            {
              if ( cur.low )
                {
                  if ((bmap.get (rx + 1, ry) || bmap.get (rx + 1, ry + 1) || bmap.get (rx, ry + 1)) &&
                      (bmap.get (rx - 1, ry - 1) || bmap.get (rx, ry - 1) || bmap.get (rx - 1, ry)))
                      bounce = true;
                }
              else
                {
                  if ((bmap.get (rx + 1, ry + 1) || bmap.get (rx, ry + 1) || bmap.get (rx + 1, ry)) &&
                      (bmap.get (rx - 1, ry) || bmap.get (rx - 1, ry - 1) || bmap.get (rx, ry - 1)))
                      bounce = true;
                }
            }
          else
            {
              if ( cur.low )
                {
                  if ((bmap.get (rx + 1, ry - 1) || bmap.get (rx, ry - 1) || bmap.get (rx + 1, ry)) &&
                      (bmap.get (rx - 1, ry) || bmap.get (rx - 1, ry + 1) || bmap.get (rx, ry + 1)))
                      bounce = true;
                }
              else
                {
                  if ((bmap.get (rx + 1, ry) || bmap.get (rx + 1, ry - 1) || bmap.get (rx, ry - 1)) &&
                      (bmap.get (rx, ry + 1) || bmap.get (rx - 1, ry + 1) || bmap.get (rx - 1, ry)))
                      bounce = true;
                }
            }

          if ( bounce )
            {
              holes++;
              // now try to detect fjords
              if ( (cur.pos.x ^ cur.pos.y) & 1 )
                {
                  if ( cur.low )
                    {
                      if ( !bmap.get (rx + 1, ry - 1) )
                          fjords++;
                    }
                  else
                    {
                      if ( !bmap.get (rx - 1, ry + 1) )
                          fjords++;
                    }
                }
              else
                {
                  if ( cur.low )
                    {
                      if ( !bmap.get (rx - 1, ry - 1) )
                          fjords++;
                    }
                  else
                    {
                      if ( !bmap.get (rx + 1, ry + 1) )
                          fjords++;
                    }
                }
            }
        }

      if ( bounce )
        {
          if ( zigzag )
            {
              if ( lc >= 6 )
                {
                  lines += lc - 5;
                  if ( corner )
                      corners++;
                  corner = true;
                }
              else
                {
                  corner = false;
                }
              lc = 1;
            }
          else
            {
              lc++;
              zigzag = true;
            }

          //bonds++;
          turns = 1;
#ifdef PRINTOUT
          //grid.push_back (cur.pos);
#endif
          cur.hamster();
#ifdef PRINTOUT
          path.push_back (cur);
#endif
        }
      else
        {
          if ( zigzag )
            {
              lc++;
              zigzag = false;
            }
          else
            {
              if ( lc >= 6 )
                {
                  lines += lc - 5;
                  if ( corner )
                      corners++;
                  corner = true;
                }
              else
                {
                  corner = false;
                }
              lc = 1;
            }

          cur.punch();
#ifdef PRINTOUT
          path.push_back (cur);
#endif
          cur.hamster();
#ifdef PRINTOUT
          path.push_back (cur);
#endif
          turns++;
          if ( turns == 4 )
              ends++;
        }
      if ( cur == end2 ) // full circle, something went wrong
          return false;
      perim++;
    }

#if 0
      coord_t spos = cur;

      hull = 0;
      bonds = 0;
      ends = 0;
      lines = 0;
      sides = 0;

      pos_t cpos, ppos, rstart, rpos, prpos;
      cpos = ppos = rstart = rpos = prpos = pos_t(0, 0);
      bool first = true;

      // now calculate the perimeter properties
      while ( 1 )
        {
//#ifndef PRINTOUT
#if 1
          if ( cur == end )
            {
              ret = true;
              break;
            }
#if 0
          if ( !first )
            {
              int x = rstart.x - rpos.x;
              int y = rstart.y - rpos.y;

              if ( abs (x) == limit / 2 || abs (y) == limit / 2 )
                {
                  ret = true;
                  break;
                }
            }
#endif
#endif

          int rx = cur.pos.x - start.pos.x;
          int ry = cur.pos.y - start.pos.y;

          bool bounce = bmap.get (rx, ry);

          if ( !bounce )
            {
              if ( (cur.pos.x ^ cur.pos.y) & 1 )
                {
                  if ( cur.low )
                    {
                      if ((bmap.get (rx + 1, ry) || bmap.get (rx + 1, ry + 1) || bmap.get (rx, ry + 1)) &&
                          (bmap.get (rx - 1, ry - 1) || bmap.get (rx, ry - 1) || bmap.get (rx - 1, ry)))
                          bounce = true;
                    }
                  else
                    {
                      if ((bmap.get (rx + 1, ry + 1) || bmap.get (rx, ry + 1) || bmap.get (rx + 1, ry)) &&
                          (bmap.get (rx - 1, ry) || bmap.get (rx - 1, ry - 1) || bmap.get (rx, ry - 1)))
                          bounce = true;
                    }
                }
              else
                {
                  if ( cur.low )
                    {
                      if ((bmap.get (rx + 1, ry - 1) || bmap.get (rx, ry - 1) || bmap.get (rx + 1, ry)) &&
                          (bmap.get (rx - 1, ry) || bmap.get (rx - 1, ry + 1) || bmap.get (rx, ry + 1)))
                          bounce = true;
                    }
                  else
                    {
                      if ((bmap.get (rx + 1, ry) || bmap.get (rx + 1, ry - 1) || bmap.get (rx, ry - 1)) &&
                          (bmap.get (rx, ry + 1) || bmap.get (rx - 1, ry + 1) || bmap.get (rx - 1, ry)))
                          bounce = true;
                    }
                }
            }

          pos_t npos = cpos;
          if ( !bounce )
            {
              if ( cur.pos.x * 2 - ppos.x == -2 && cur.pos.y * 2 == ppos.y )
                  npos = pos_t(cur.pos.x * 2 + 1, cur.pos.y * 2 - 1);
              else if ( cur.pos.x * 2 - ppos.x == 2 && cur.pos.y * 2 == ppos.y )
                  npos = pos_t(cur.pos.x * 2 - 1, cur.pos.y * 2 + 1);
              else if ( cur.pos.y * 2 - ppos.y == -2 && cur.pos.x * 2 == ppos.x )
                  npos = pos_t(cur.pos.x * 2 + 1, cur.pos.y * 2 + 1);
              else if ( cur.pos.y * 2 - ppos.y == 2 && cur.pos.x * 2 == ppos.x )
                  npos = pos_t(cur.pos.x * 2 - 1, cur.pos.y * 2 - 1);
              ppos = pos_t(cur.pos.x * 2, cur.pos.y * 2);
            }
          else
            {
              if ( (cur.pos.x ^ cur.pos.y) & 1 )
                {
                  if ( cur.low )
                      npos = pos_t(cur.pos.x * 2 - 1, cur.pos.y * 2 + 1);
                  else
                      npos = pos_t(cur.pos.x * 2 + 1, cur.pos.y * 2 - 1);
                }
              else
                {
                  if ( cur.low )
                      npos = pos_t(cur.pos.x * 2 + 1, cur.pos.y * 2 + 1);
                  else
                      npos = pos_t(cur.pos.x * 2 - 1, cur.pos.y * 2 - 1);
                }
              ppos = npos;
            }

          if ( cpos != npos )
            {
              if ( bounce )
                  sides++;
              pos_t nrpos = pos_t((npos.x + npos.y) / 4 - start.pos.x, (npos.y - npos.x) / 4);
              if ( hull >= 1 && ((nrpos.x == rpos.x && prpos.x == rpos.x) || (nrpos.y == rpos.y && prpos.y == rpos.y)) )
                  lines++;
              prpos = rpos;
              rpos = nrpos;
              if ( first )
                {
                  rstart = rpos;
                  first = false;
                }
              else
                {
                  hull++;
                }
              cpos = npos;
#ifdef PRINTOUT
              perimpath.push_back (npos);
#endif
            }
              

          if ( bounce )
            {
              cur.hamster();
            }
          else
            {
              cur.punch();
              cur.hamster();
            }
          if ( cur == spos )
              break;
        }
#endif

#ifdef PRINTOUT
  printf ("Printout: printing now, having taken %zu hull steps\n", hull);
  print_path ("printout_path.csv", path);
  print_path ("printout_path1.csv", path1);
  print_path ("printout_path2.csv", path2);
  print_grid ("printout_grid.csv", grid);
  //print_perimpath ("printout_perimpath.csv", perimpath);
  print_surf (grid);
  exit(0);
#endif

  return ret;
}

bool Walker::walk_perim (int limit)
{
  //int limit2 = limit * limit;
  bool ret = false;
  coord_t cur = start;

  static Bitmap bmap(32);
  int search_limit = limit * 2;
  bmap.reset (search_limit + 1);

#ifdef PRINTOUT
  path_t path;
  std::vector<pos_t> perimpath, mypath;
  grid_t grid;
  path.push_back (cur);
#endif

  // at first find whole clusters of desired size
  size_t hs = 0;
  int rc = 0;
  int maxx = 0, minx = FIELD_SIZE, maxy = 0, miny = FIELD_SIZE;
  while ( 1 )
    {
      int ox = start.pos.x - cur.pos.x;
      int oy = start.pos.y - cur.pos.y;
      int x = ox;
      int y = oy;

      maxx = std::max (maxx, cur.pos.x);
      maxy = std::max (maxy, cur.pos.y);
      minx = std::min (minx, cur.pos.x);
      miny = std::min (miny, cur.pos.y);

      if ( abs (x) == search_limit || abs (y) == search_limit )
          return false;

      if ( field->valueb (cur.pos.x, cur.pos.y) )
        {
#ifdef PRINTOUT
          grid.push_back (cur.pos);
#endif
          bmap.set (cur.pos.x - start.pos.x, cur.pos.y - start.pos.y);
          cur.hamster();
#ifdef PRINTOUT
          path.push_back (cur);
#endif
          rc--;
        }
      else
        {
          cur.punch();
#ifdef PRINTOUT
          path.push_back (cur);
#endif
          cur.hamster();
#ifdef PRINTOUT
          path.push_back (cur);
#endif
          rc++;
        }
      if ( cur == start )
          break;
      hs++;
    }

  // inside perimeter
  if ( rc < 0 || ((maxx - minx) < limit && (maxy - miny) < limit) )
      return false;

  //printf ("Found it, hs=%zu, rc=%d!\n", hs, rc);

  size_t target = 0;
  size_t curhs = 0;
  coord_t spos = start;

  // try different position on the hull, find one that is not boxed and of sufficient size
  while ( target < hs ) 
    {
      // skip on some of the way to find a better starting position
      cur = spos;
      while ( curhs < target )
        {
          int rx = cur.pos.x - start.pos.x;
          int ry = cur.pos.y - start.pos.y;

          if ( bmap.get (rx, ry) )
            {
              cur.hamster();
            }
          else
            {
              cur.punch();
              cur.hamster();
            }

          curhs++;
          if ( curhs > hs ) return false;
        }
      spos = cur;

      hull = 0;
      bonds = 0;
      ends = 0;
      lines = 0;
      sides = 0;

      pos_t cpos, ppos, rstart, rpos, prpos;
      cpos = ppos = rstart = rpos = prpos = pos_t(0, 0);
      bool first = true;

      // now calculate the perimeter properties
      while ( 1 )
        {
#ifndef PRINTOUT
          if ( !first )
            {
              int x = rstart.x - rpos.x;
              int y = rstart.y - rpos.y;

              if ( abs (x) == limit || abs (y) == limit )
                {
                  ret = true;
                  break;
                }
            }
#endif

          int rx = cur.pos.x - start.pos.x;
          int ry = cur.pos.y - start.pos.y;

          bool bounce = bmap.get (rx, ry);

          if ( !bounce )
            {
              if ( (cur.pos.x ^ cur.pos.y) & 1 )
                {
                  if ( cur.low )
                    {
                      if ((bmap.get (rx + 1, ry) || bmap.get (rx + 1, ry + 1) || bmap.get (rx, ry + 1)) &&
                          (bmap.get (rx - 1, ry - 1) || bmap.get (rx, ry - 1) || bmap.get (rx - 1, ry)))
                          bounce = true;
                    }
                  else
                    {
                      if ((bmap.get (rx + 1, ry + 1) || bmap.get (rx, ry + 1) || bmap.get (rx + 1, ry)) &&
                          (bmap.get (rx - 1, ry) || bmap.get (rx - 1, ry - 1) || bmap.get (rx, ry - 1)))
                          bounce = true;
                    }
                }
              else
                {
                  if ( cur.low )
                    {
                      if ((bmap.get (rx + 1, ry - 1) || bmap.get (rx, ry - 1) || bmap.get (rx + 1, ry)) &&
                          (bmap.get (rx - 1, ry) || bmap.get (rx - 1, ry + 1) || bmap.get (rx, ry + 1)))
                          bounce = true;
                    }
                  else
                    {
                      if ((bmap.get (rx + 1, ry) || bmap.get (rx + 1, ry - 1) || bmap.get (rx, ry - 1)) &&
                          (bmap.get (rx, ry + 1) || bmap.get (rx - 1, ry + 1) || bmap.get (rx - 1, ry)))
                          bounce = true;
                    }
                }
            }

          pos_t npos = cpos;
          if ( !bounce )
            {
              if ( cur.pos.x * 2 - ppos.x == -2 && cur.pos.y * 2 == ppos.y )
                  npos = pos_t(cur.pos.x * 2 + 1, cur.pos.y * 2 - 1);
              else if ( cur.pos.x * 2 - ppos.x == 2 && cur.pos.y * 2 == ppos.y )
                  npos = pos_t(cur.pos.x * 2 - 1, cur.pos.y * 2 + 1);
              else if ( cur.pos.y * 2 - ppos.y == -2 && cur.pos.x * 2 == ppos.x )
                  npos = pos_t(cur.pos.x * 2 + 1, cur.pos.y * 2 + 1);
              else if ( cur.pos.y * 2 - ppos.y == 2 && cur.pos.x * 2 == ppos.x )
                  npos = pos_t(cur.pos.x * 2 - 1, cur.pos.y * 2 - 1);
              ppos = pos_t(cur.pos.x * 2, cur.pos.y * 2);
            }
          else
            {
              if ( (cur.pos.x ^ cur.pos.y) & 1 )
                {
                  if ( cur.low )
                      npos = pos_t(cur.pos.x * 2 - 1, cur.pos.y * 2 + 1);
                  else
                      npos = pos_t(cur.pos.x * 2 + 1, cur.pos.y * 2 - 1);
                }
              else
                {
                  if ( cur.low )
                      npos = pos_t(cur.pos.x * 2 + 1, cur.pos.y * 2 + 1);
                  else
                      npos = pos_t(cur.pos.x * 2 - 1, cur.pos.y * 2 - 1);
                }
              ppos = npos;
            }

          if ( cpos != npos )
            {
              if ( bounce )
                  sides++;
              pos_t nrpos = pos_t((npos.x + npos.y) / 4 - start.pos.x, (npos.y - npos.x) / 4);
              if ( hull >= 1 && ((nrpos.x == rpos.x && prpos.x == rpos.x) || (nrpos.y == rpos.y && prpos.y == rpos.y)) )
                  lines++;
              prpos = rpos;
              rpos = nrpos;
              if ( first )
                {
                  rstart = rpos;
                  first = false;
                }
              else
                {
                  hull++;
                }
              cpos = npos;
#ifdef PRINTOUT
              //mypath.push_back (npos);
              perimpath.push_back (npos);
              mypath.push_back (pos_t(rpos.x * 2, rpos.y * 2));
#endif
            }
              

          if ( bounce )
            {
              cur.hamster();
            }
          else
            {
              cur.punch();
              cur.hamster();
            }
          if ( cur == spos )
              break;
        }

      if ( ret )
          return true;

      target += hs / 8;
#ifdef PRINTOUT
      if ( !path.empty() )
          break;
#endif
    }

#ifdef PRINTOUT
  printf ("Printout: printing now, having taken %zu hull steps\n", hull);
  print_path ("printout_path.csv", path);
  print_grid ("printout_grid.csv", grid);
  print_perimpath ("printout_perimpath.csv", perimpath);
  print_perimpath ("printout_mypath.csv", mypath);
  print_surf (grid);
  exit(0);
#endif
  return ret;
}

#define SPACERATIO 0.2

void Walker::print_path (const char *filen, const path_t& path)
{
  assert (path.size() > 1);
  FILE *fp = fopen (filen, "w");
  for ( size_t i = 0; i < path.size(); i++ )
    {
      coord_t cur = path[i];
      vect2d pos(cur.pos.x, cur.pos.y);

      if ( i == 0 )
        {
          vect2d npos(path[i + 1].pos.x, path[i + 1].pos.y);
          vect2d p = pos + (npos - pos) * SPACERATIO;
          fprintf (fp, "%.2f %.2f\n", p.x, p.y);
          continue;
        }

      vect2d ppos(path[i - 1].pos.x, path[i - 1].pos.y);

      if ( path[i - 1].pos != cur.pos )
        {
          vect2d dir = pos - ppos;
          vect2d p = ppos + dir * SPACERATIO;
          fprintf (fp, "%.2f %.2f\n", p.x, p.y);
          p = pos - dir * SPACERATIO;
          fprintf (fp, "%.2f %.2f\n", p.x, p.y);
        }
#if 0
      vect2d box;
      if ( cur.pos.x & 1 )
        {
          if ( cur.pos.y & 1 )
            {
              if ( cur.low )
                  box = vect2d (pos.x + 0.5, pos.y + 0.5);
              else
                  box = vect2d (pos.x - 0.5, pos.y - 0.5);
            }
          else
            {
              if ( cur.low )
                  box = vect2d (pos.x - 0.5, pos.y + 0.5);
              else
                  box = vect2d (pos.x + 0.5, pos.y - 0.5);
            }
        }
      else
        {
          if ( cur.pos.y & 1 )
            {
              if ( cur.low )
                  box = vect2d (pos.x - 0.5, pos.y + 0.5);
              else
                  box = vect2d (pos.x + 0.5, pos.y - 0.5);
            }
          else
            {
              if ( cur.low )
                  box = vect2d (pos.x + 0.5, pos.y + 0.5);
              else
                  box = vect2d (pos.x - 0.5, pos.y - 0.5);
            }
        }
      vect2d dir = (box - pos).norm();
      vect2d rpos = pos + dir * 0.1;
      fprintf (fp, "%.2f %.2f\n", rpos.x, rpos.y);
#endif
    }
  fclose (fp);
}

void Walker::print_perimpath (const char *filen, const std::vector<pos_t>& path)
{
  assert (!path.empty());
  FILE *fp = fopen (filen, "w");
  for ( size_t i = 0; i < path.size(); i++ )
    {
      fprintf (fp, "%f %f\n", (double)path[i].x / 2, (double)path[i].y / 2);
    }
  //fprintf (fp, "%f %f\n", (double)path[0].x / 2, (double)path[0].y / 2);
  fclose (fp);
}

#ifdef PRINTOUT
  #if __cplusplus >= 201103L
    #include <unordered_set>
  #endif
#endif

void Walker::print_grid (const char *path, const grid_t& grid)
{
#ifdef PRINTOUT
#if __cplusplus >= 201103L
  std::unordered_set<pos_t, pos_t::hash_t> dups;
  FILE *fp = fopen (path, "w");
  for ( size_t i = 0; i < grid.size(); i++ )
    {
      pos_t pos = grid[i];
      if ( dups.find (pos) != dups.end() )
          continue;
      dups.insert (pos);
      vect2d vpos (pos.x, pos.y);
      vect2d dir;
      if ( (pos.x ^ pos.y) & 1 )
          dir = vect2d (0.5, 0.5);
      else
          dir = vect2d (0.5, -0.5);

      if ( i == 0 ) // hack to throw off isolines
          fprintf (fp, "%.2f %.2f\n", (vpos + dir).x, (vpos + dir).y);
      fprintf (fp, "%.2f %.2f\n", (vpos + dir).x, (vpos + dir).y);
      fprintf (fp, "%.2f %.2f\n", (vpos - dir).x, (vpos - dir).y);
      fprintf (fp, "\n");
    }
  fclose (fp);
#else
  fprintf(stderr, "Walker::print_grid: not supported, need c++11!\n");
#endif
#endif
}

void Walker::print_surf (const grid_t& grid)
{
  assert (grid.size() > 0);
  pos_t ll, ur;
  ll = ur = grid[0];
  for ( size_t i = 0; i < grid.size(); i++ )
    {
      if ( grid[i].x < ll.x )
          ll.x = grid[i].x;
      if ( grid[i].x > ur.x )
          ur.x = grid[i].x;
      if ( grid[i].y < ll.y )
          ll.y = grid[i].y;
      if ( grid[i].y > ur.y )
          ur.y = grid[i].y;
    }
  ll.x -= 8;
  ll.y -= 8;
  ur.x += 8;
  ur.y += 8;

  FILE *fp = fopen ("printout_surf.csv", "w");
  for ( int x = ll.x; x <= ur.x; x++ )
    {
      for ( int y = ll.y; y <= ur.y; y++ )
        {
          double v = v2d(field->value (x, y));
          //v = v >= 0 ? 0 : 1;
          v = -v;
          fprintf (fp, "%d %d %f\n", x, y, v);
        }
      fprintf (fp, "\n");
    }
  fclose (fp);
}

void freq::print (const char *basename, int A)
{
  const char *dp = strrchr (basename, '.');
  if ( !dp )
      dp = basename + strlen (basename);
  std::string basepath(basename, dp);
  basepath += ".";
  std::ostringstream oss;
  oss << A;
  basepath += oss.str();
  basepath += ".csv";

  FILE *fp = fopen (basepath.c_str(), "w");
  for ( f_t::iterator it = f.begin(); it != f.end(); it++ )
    {
      fprintf (fp, "%d %d\n", it->first, it->second);
    }
  fclose (fp);
}

#if 0
#define DIST_DIV 200
#define DIST_RANGE 4.0
#define DIST_COUNT 1000
#define DIST_SIDE ((int)(DIST_DIV*DIST_RANGE))
#define DIST_START_LEVEL 17
void Field::calc_dist (const char *path, size_t samples, size_t levels)
{
  double fact[levels + 1];
  for ( size_t j = 0; j <= levels; j++ )
      fact[j] = sqrtl (ALPHA * ((long double)1.0 - powl (BASEPOW, -EXPONENT))) * powl (BASEPOW, -EXPONENT * ((long double)j / 2.0));

  size_t count[2 * DIST_SIDE + 1];
  memset (count, 0, sizeof(count));

  long double totcum2 = 0;
  size_t ncum = 0;

  for ( size_t s = 0; s < samples; s++ )
    {
      for ( size_t i = 0; i < DIST_COUNT; i++ )
        {
          long double cum = 0;
          for ( size_t j = DIST_START_LEVEL; j <= levels; j++ )
              cum += calc_cell() * fact[j];
          totcum2 += cum * cum;
          ncum++;
          int index = round(cum * DIST_DIV);
          if ( index < -DIST_SIDE || index > DIST_SIDE )
              continue;
          count[DIST_SIDE + index]++;
        }
      fprintf (stderr, "flush!\n");
      //break;
    }

  FILE *fp = fopen (path, "w");
  for ( int i = -DIST_SIDE; i <= DIST_SIDE; i++ )
    {
      if ( !count[i + DIST_SIDE] )
          continue;
      fprintf (fp, "%g %zu\n", (double)i / DIST_DIV, count[i + DIST_SIDE]);
    }
  fclose (fp);
  {
    const char *basename = strrchr (path, '/');
    if ( basename )
        basename++;
    else
        basename = path;
    std::string newpath = std::string(path, path + (basename - path)) + "var." + basename;
    FILE *fp = fopen (newpath.c_str(), "w");
    fprintf (fp, "1 %zu\n", ncum);
    fprintf (fp, "2 %.12Le\n", totcum2);
    fclose (fp);
  }
}

#undef DIST_COUNT
#define DIST_COUNT 1000000000

void Field::calc_cell_dist (const char *path)
{
  size_t count[2 * DIST_SIDE + 1];
  memset (count, 0, sizeof(count));

  long double totalsum = 0;
  long double totalsum2 = 0;

  for ( size_t i = 0; i < DIST_COUNT; i++ )
    {
      double val = calc_cell();
      totalsum += val;
      totalsum2 += val * val;
      int index = round(val * DIST_DIV);
      if ( index < -DIST_SIDE || index > DIST_SIDE )
          continue;
      count[DIST_SIDE + index]++;
    }
  fprintf (stderr, "E(X)=%g, var(X)=%g\n", (double)(totalsum / DIST_COUNT), (double)(totalsum2 / DIST_COUNT));
  FILE *fp = fopen (path, "w");
  for ( int i = -DIST_SIDE; i <= DIST_SIDE; i++ )
    {
      if ( !count[i + DIST_SIDE] )
          continue;
      fprintf (fp, "%g %zu\n", (double)i / DIST_DIV, count[i + DIST_SIDE]);
    }
  fclose (fp);

  {
  const char *basename = strrchr (path, '/');
  if ( basename )
      basename++;
  else
      basename = path;
  std::string newpath = std::string(path, path + (basename - path)) + "var." + basename;
  FILE *fp = fopen (newpath.c_str(), "w");
  fprintf (fp, "1 %zu\n", (size_t)DIST_COUNT);
  fprintf (fp, "2 %.12Le\n", totalsum2);
  fclose (fp);
  }
}
#endif

#ifndef RNDF_MMAP
RNDF::RNDF ()
{
  rng = gsl_rng_alloc (gsl_rng_default);
  fp = 0;
  seed (gsl_rng_default, 0xffffffffffffffffL);
}

RNDF::~RNDF ()
{
  gsl_rng_free (rng);
  if ( fp )
      fclose (fp);
}

void RNDF::seed (const gsl_rng_type *type, unsigned long int n, const char *binpath)
{
  if ( fp )
      fclose (fp);
  if ( !binpath )
      binpath = "rnd.bin";
  fp = fopen (binpath, "rb");
  if ( !fp )
    {
      printf ("Could not find random data from '%s', falling back to pseudorandom numbers!\n", binpath);
    }
  else
    {
      struct stat st;
      assert (stat (binpath, &st) == 0);
      flen = st.st_size;
      assert (flen < 0x7fffffff);
      assert (flen > 10 * RNDF_BLOCK_SIZE);
      printf ("Opened binary data '%s', file length = %zu.\n", binpath, flen);
    }

  printf ("Using random generator: %s\n", type->name);
  if ( !n )
    {
      printf ("Seed given as 0, generating from %s\n", RANDOM_DEVICE);
      FILE *rndin = fopen (RANDOM_DEVICE, "rb");
      if ( !rndin )
        {
          printf ("Failed to open %s!\n", RANDOM_DEVICE);
          exit(-1);
        }
      fread (&n, sizeof(n), 1, rndin);
      fclose (rndin);
    }
  printf ("Using seed %lu\n", n);
  gsl_rng_free (rng);
  rng = gsl_rng_alloc (type);
  gsl_rng_set (rng, n);
  for ( int i = 0; i < 8; i++ )
      b.bits[i] = gsl_rng_uniform_int (rng, 256);
  int exponent = 1023;
#ifdef BIGENDIAN
  b.bits[0] = exponent >> 4;
  b.bits[1] = (b.bits[1] & 0x0f) | ((exponent & 0x0f) << 4);
#else
  b.bits[7] = exponent >> 4;
  b.bits[6] = (b.bits[6] & 0x0f) | ((exponent & 0x0f) << 4);
#endif
  reload();
}

double RNDF::next ()
{
  if ( !fp )
      return gsl_rng_uniform (rng);

  if ( dptr - data == RNDF_BLOCK_SIZE )
      reload();
#ifdef BIGENDIAN
  b.bits[1] = (b.bits[1] & (0xf0)) | ((b.bits[6] ^ *dptr++) & 0x0f);
  b.bits[2] ^= *dptr++;
  b.bits[3] ^= *dptr++;
  b.bits[4] ^= *dptr++;
  b.bits[5] ^= *dptr++;
  b.bits[6] ^= *dptr++;
  b.bits[7] ^= *dptr++;
#else
  b.bits[0] ^= *dptr++;
  b.bits[1] ^= *dptr++;
  b.bits[2] ^= *dptr++;
  b.bits[3] ^= *dptr++;
  b.bits[4] ^= *dptr++;
  b.bits[5] ^= *dptr++;
  b.bits[6] = (b.bits[6] & (0xf0)) | ((b.bits[6] ^ *dptr++) & 0x0f);
#endif
  return b.dbl - 1.0;
}

void RNDF::reload ()
{
  if ( !fp ) return;
  int offset1, offset2;

  while (1)
    {
      offset1 = gsl_rng_uniform_int (rng, flen - RNDF_BLOCK_SIZE);
      offset2 = gsl_rng_uniform_int (rng, flen - RNDF_BLOCK_SIZE);
      if ( abs(offset1 - offset2) >= RNDF_BLOCK_SIZE )
          break;
    }

  fseek (fp, offset1, SEEK_SET);
  assert (fread (block1, RNDF_BLOCK_SIZE, 1, fp) == 1);
  fseek (fp, offset2, SEEK_SET);
  assert (fread (block2, RNDF_BLOCK_SIZE, 1, fp) == 1);
  for ( size_t i = 0; i < RNDF_BLOCK_SIZE; i++ )
      data[i] = block1[i] ^ block2[i];
  dptr = data;
}

#else

RNDF::RNDF ()
{
  rng = gsl_rng_alloc (gsl_rng_default);
  memfd = -1;
  seed (gsl_rng_default, 0xffffffffffffffffL);
}

RNDF::~RNDF ()
{
  gsl_rng_free (rng);
  if ( memfd != -1 )
    {
      munmap (mem, flen);
      close (memfd);
    }
}

void RNDF::seed (const gsl_rng_type *type, unsigned long int n, const char *binpath)
{
  if ( memfd != -1 )
    {
      munmap (mem, flen);
      close (memfd);
    }
  if ( !binpath )
      binpath = "rnd.bin";
  memfd = open (binpath, O_RDONLY);
  if ( memfd == -1 )
    {
      printf ("Could not find random data from '%s', falling back to pseudorandom numbers!\n", binpath);
    }
  else
    {
      struct stat st;
      assert (stat (binpath, &st) == 0);
      flen = st.st_size;
      assert (flen < 0x7fffffff);
      assert (flen > 10 * RNDF_BLOCK_SIZE);
      printf ("Opened binary data '%s', file length = %zu.\n", binpath, flen);
      mem = (unsigned char *)mmap (NULL, st.st_size, PROT_READ, MAP_SHARED|MAP_POPULATE, memfd, 0);
      //mem = (unsigned char *)mmap (NULL, st.st_size, PROT_READ, MAP_SHARED|MAP_POPULATE, memfd, 0);
      assert (mem);
    }

  printf ("Using random generator: %s\n", type->name);
  if ( !n )
    {
      printf ("Seed given as 0, generating from %s\n", RANDOM_DEVICE);
      FILE *rndin = fopen (RANDOM_DEVICE, "rb");
      if ( !rndin )
        {
          printf ("Failed to open %s!\n", RANDOM_DEVICE);
          exit(-1);
        }
      fread (&n, sizeof(n), 1, rndin);
      fclose (rndin);
    }
  printf ("Using seed %lu\n", n);
  gsl_rng_free (rng);
  rng = gsl_rng_alloc (type);
  gsl_rng_set (rng, n);
  for ( int i = 0; i < 8; i++ )
      b.bits[i] = gsl_rng_uniform_int (rng, 256);
  int exponent = 1023;
#ifdef BIGENDIAN
  b.bits[0] = exponent >> 4;
  b.bits[1] = (b.bits[1] & 0x0f) | ((exponent & 0x0f) << 4);
#else
  b.bits[7] = exponent >> 4;
  b.bits[6] = (b.bits[6] & 0x0f) | ((exponent & 0x0f) << 4);
#endif
  reload();
}

double RNDF::next ()
{
  if ( memfd == -1 )
      return gsl_rng_uniform (rng);

  if ( dptr - data == RNDF_BLOCK_SIZE )
      reload();
#ifdef BIGENDIAN
  b.bits[1] = (b.bits[1] & (0xf0)) | ((b.bits[6] ^ *dptr++) & 0x0f);
  b.bits[2] ^= *dptr++;
  b.bits[3] ^= *dptr++;
  b.bits[4] ^= *dptr++;
  b.bits[5] ^= *dptr++;
  b.bits[6] ^= *dptr++;
  b.bits[7] ^= *dptr++;
#else
  b.bits[0] ^= *dptr++;
  b.bits[1] ^= *dptr++;
  b.bits[2] ^= *dptr++;
  b.bits[3] ^= *dptr++;
  b.bits[4] ^= *dptr++;
  b.bits[5] ^= *dptr++;
  b.bits[6] = (b.bits[6] & (0xf0)) | ((b.bits[6] ^ *dptr++) & 0x0f);
#endif
  return b.dbl - 1.0;
}

void RNDF::reload ()
{
  if ( memfd == -1 ) return;
  int offset1, offset2;

  while (1)
    {
      offset1 = gsl_rng_uniform_int (rng, flen - RNDF_BLOCK_SIZE);
      offset2 = gsl_rng_uniform_int (rng, flen - RNDF_BLOCK_SIZE);
      if ( abs(offset1 - offset2) >= RNDF_BLOCK_SIZE )
          break;
    }

  unsigned char *b1 = mem + offset1;
  unsigned char *b2 = mem + offset2;
  for ( size_t i = 0; i < RNDF_BLOCK_SIZE; i++ )
      data[i] = *b1++ ^ *b2++;
  dptr = data;
}
#endif

